# Source code for BigDFT.UnitCells

```"""
A module for handling unit cells for periodic calculations.
"""

[docs]class UnitCell:
"""
Defines a wrapper for unit cells.

Args:
cell (list): a list of unit cell vectors. This can either be a 3x1 list
or a 3x3 list. Currently, only orthorhombic cells are supported.
units (str): the units of the cell parameters. If cell is set to
None, the free boundary condition is enforced.
units (str): the unit of length.
"""
def __init__(self, cell=None, units="bohr"):
from BigDFT.Atoms import AU_to_A, IsAngstroem

if cell is None:
self.cell = [[float("inf"), 0, 0],
[0, float("inf"), 0],
[0, 0, float("inf")]]
return

# Copy cell based on the format
if isinstance(cell, list):
self.cell = cell
else:
self.cell = [[cell, 0, 0], [0, cell, 0], [0, 0, cell]]

# Check that the unit cell is valid
if self[2, 2] == float("inf") and not (self[1, 1] == float("inf") and
self[0, 0] == float("inf")):
raise ValueError("Boundary condition: ",
"If z is infinity, x and y must also be.")
elif self[0, 0] == float("inf") and self[1, 1] != float("inf"):
raise ValueError("Boundary condition: ",
"If x is infinity, y must also be.")
elif self[0, 1] != 0 or self[0, 2] != 0 or \
self[1, 0] != 0 or self[1, 2] != 0 or \
self[2, 0] != 0 or self[2, 1] != 0:
raise ValueError("Non-orthorhombic cells not supported")

for i in range(3):
for j in range(3):
if self[i, j] < 0:
raise ValueError("Unit cell must be non-negative values.")

# Store internally in bohr
if IsAngstroem(units):
for i in range(3):
for j in range(3):
self.cell[i][j] /= AU_to_A

def __getitem__(self, idx):
return self.cell[idx][idx]

[docs]    def get_boundary_condition(self, units="bohr"):
"""
Get a string description of the boundary condition (i.e. free,
surface, wire, periodic)

Args:
units (str): the units to report the cell in.

Returns:
(str): a string description of the boundary condition.
"""
from BigDFT.Atoms import AU_to_A, IsAngstroem

# Match units
if IsAngstroem(units):
conv = AU_to_A
else:
conv = 1

if all([self[i, i] == float("inf") for i in range(3)]):
return "free"
elif self[0, 0] == float("inf") and self[1, 1] == float("inf"):
return "wire 0.0 0.0 " + str(self[2, 2]*conv)
elif self[1, 1] == float("inf"):
return "surface " + str(self[0, 0]*conv) + " 0.0 " + \
str(self[2, 2]*conv)
else:
return "periodic " + " ".join([str(conv*self[i, i])
for i in range(3)])

[docs]    def get_posinp(self, units="bohr"):
"""
Create the dictionary representation of the cell that is passed to
BigDFT.

Returns:
(list): a list of the three values of the unit cell.
"""
from BigDFT.Atoms import AU_to_A, IsAngstroem, IsBohr
if IsAngstroem(units):
return [self[i, i] * AU_to_A for i in range(3)]
elif IsBohr(units):
return [self[i, i] for i in range(3)]
else:
raise ValueError("Posinp units must be bohr or angstroem")

[docs]    def minimum_image(self, pos, units="bohr"):
"""
Given a vector of three positions, this wraps those positions inside
the cell using the minimum image convention.

Returns:
(list): a list of the values of the wrapped position.
"""
from BigDFT.Atoms import AU_to_A, IsAngstroem, IsReduced, IsBohr

# Match units
if IsAngstroem(units):
conversion = 1/AU_to_A
elif IsReduced(units):
return pos
elif IsBohr(units):
conversion = 1
else:
raise ValueError("Invalid unit: " + units)

bohrpos = [x*conversion for x in pos]

for i in range(3):
if self[i, i] == float("inf"):
continue
while(bohrpos[i] > self[i, i]):
bohrpos[i] -= self[i, i]
while(bohrpos[i] < 0):
bohrpos[i] += self[i, i]

# Convert back
return [x/conversion for x in bohrpos]

[docs]    def to_cartesian(self, pos):
"""
Convert a vector which is in reduced units to cartesian units.

Returns:
(list): the position in cartesian coordinates.
"""
from numpy import array
return array(self.cell).dot(pos)

[docs]    def to_reduced(self, pos):
"""
Convert a vector which is in reduced units to cartesian units.

Returns:
(list): the position in reduced coordinates.
"""
from numpy import array
from numpy.linalg import solve

return solve(array(self.cell), pos)

def _example():
# Create a basic unit cell
cell = UnitCell([10, 8, 4], units="angstroem")

# Print out the posinp representation
print(cell.get_posinp())
print(cell.get_posinp(units="angstroem"))

# Right now we enforce the orthorhombic condition
try:
cell = UnitCell([[10, 0, 0], [0, 8, 0], [0, 4, 10]])
except ValueError as e:
print(e)
cell = UnitCell([[10, 0, 0], [0, 8, 0], [0, 0, 4]])
print(cell.get_boundary_condition("angstroem"))

# Wire boundary condition
wire = UnitCell([float("inf"), float("inf"), 4])
print(wire.get_posinp())
print(wire.get_boundary_condition())

# Surface boundary condition
surface = UnitCell([10, float("inf"), 4])
print(surface.get_posinp())
print(surface.get_boundary_condition())

# Wrap positions to the minimum image convention.
pos = [-5, -2, -3]
print(cell.minimum_image(pos))
print(wire.minimum_image(pos))
print(surface.minimum_image(pos))

pos = [15, 12, 13]
print(cell.minimum_image(pos))
print(wire.minimum_image(pos))
print(surface.minimum_image(pos))

if __name__ == "__main__":
_example()
```